Reviewer 1: Given the extensive collection of metabolomics data (targeted and untargeted) for so many tissues and temporal timepoints (Figure 1C) it will be interesting to explore the changes of several key cellular metabolites in addition to the KEGG pathways. For example, it will be interesting to see sex- and time-specific changes in the plasma, muscles, heart of such metabolites as glucose, pyruvate, lactate, acetate, perhaps key intermediates of glycolysis and the TCA cycle. If the data is available, it will be also interesting to characterize the changes in the energy and redox ratios, i.e. ATP/ADP, NADH/NAD+ (or individual concentrations, such as ATP). Given the central roles of these metabolites in cellular physiology/health and multiple changes in the corresponding metabolic pathways (Figure 7), it will be good to present and discuss the observed metabolic changes.

library(MotrpacRatTraining6mo) # v1.5.2
library(data.table)
library(ggplot2)
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
knitr::opts_knit$set(root.dir = "~/Documents/motrpac/MAWG/MANUSCRIPT/for_revisions/differential_metabolites/")

Read in table of training-regulated metabolites we want to plot. From: https://docs.google.com/spreadsheets/d/1OuZWnUdg1GkTHmGFFcar01t8MQy7a0N4EAOxSPOs4sA/edit?usp=sharing

# from 
v1 = fread("differential_metabolites_1.txt")
v2 = fread("differential_metabolites_2.txt")

First, let’s plot individual metabolites. Each of these metabolites are significantly training-regulated at 5% FDR. These were manually selected as potentially interesting metabolites to highlight from the full list of differential metabolites. For each metabolite, the trajectory of the normalized sample-level data is shown on the left, and the log fold-changes from the timewise differential analysis are shown on the right. Multiple traces are shown if multiple platforms measured the metabolite.

metabolites = v1[plot=="X", feature]

for(m in metabolites){
    
  # sample-level data
  p = plot_feature_normalized_data(feature=m, add_adj_p=TRUE, facet_by_sex=TRUE)
  if(!is.null(p)){
    print(p)
  }else{
    next
  }
  
  # log fold-changes
  p = plot_feature_logfc(feature=m, facet_by_sex=TRUE)
  if(!is.null(p)){
    print(p)
  }
}

Now, let’s plot some metabolites that are differential in many tissues. We include metabolites from above that were training-regulated in multiple tissues as well as all metabolites that are training-regulated in at least 8 tissues. Trajectories of log fold-changes from the timewise differential analysis are shown, with traces colored by tissue. The training adjusted p-values, which were used to select differential features at 5% FDR, are indicated in the legend.

featureIDs = unique(c(v1[plot=="X", feature_ID], v2[plot=="X", feature_ID]))
for(f in featureIDs){
  # get all training-regulated features
  features = unique(TRAINING_REGULATED_FEATURES$feature[TRAINING_REGULATED_FEATURES$feature_ID==f])
  if(length(features) <= 1){
    next
  }
  # get all logFCs
  data = list()
  for(feat in features){
    data[[feat]] = plot_feature_logfc(feature=feat, add_adj_p = TRUE, return_data = TRUE)
  }
  data = rbindlist(data)
  data[,plotting_group := paste0(tissue,feature)]
  # add adj p to legend
  data[,tissue_p := sprintf("%s (%s)", tissue, signif(selection_fdr, digits=3))]
  labels = unique(data[,tissue_p])
  names(labels) = gsub(" .*","",labels)
  
  data_df = data.frame(data)
  g = ggplot2::ggplot(data_df, ggplot2::aes(y=logFC, x=comparison_group, group=plotting_group, color=tissue)) +
    ggplot2::geom_point(size=2, position=ggplot2::position_dodge(width=0.3)) +
    ggplot2::geom_line(position=ggplot2::position_dodge(width=0.3)) +
    ggplot2::geom_errorbar(ggplot2::aes(ymin=logFC-logFC_se, ymax=logFC+logFC_se),
                           width=0.2, 
                           position=ggplot2::position_dodge(width=0.3)) +
    ggplot2::theme_classic() +
    ggplot2::geom_hline(yintercept = 0,linetype="dotted") +
    ggplot2::facet_wrap(~sex) +
    ggplot2::labs(title=f, x="Time trained (weeks)", y="Log fold-change") +
    ggplot2::theme(plot.title = ggplot2::element_text(hjust=0.5),
                   plot.subtitle = ggplot2::element_text(hjust=0.5),
                   panel.grid.major = ggplot2::element_blank(),
                   panel.grid.minor = ggplot2::element_blank()) +
    ggplot2::scale_colour_manual(values=TISSUE_COLORS[names(TISSUE_COLORS) %in% data[,tissue]], 
                                 name="Tissue",
                                 labels=labels) +
    scale_x_discrete(limits=c('control','1w','2w','fill','4w',rep('fill',3), '8w'),
                     labels=c('0','1','2','4','8'),
                     breaks=c('control','1w','2w','4w','8w'))
  print(g)
}